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ABSTRACT 

We propose a simple analytic model for the innermost (within the light 
cylinder of canonical radius, ~ c/Vt) structure of open-magnetic-field lines of 
a rotating neutron star (NS) with relativistic outflow of charged particles (elec- 
trons/positrons) and arbitrary angle between the NS spin and magnetic axes. 
We present the self-consistent solution of Maxwell’s equations for the magnetic 
field and electric current in the pair-starved regime where the density of electron- 
positron plasma generated above the pulsar polar cap is not sufficient to com- 
pletely screen the accelerating electric field and thus establish the E • B = 0 
condition above the pair-formation front up to the very high altitudes within 
the light cylinder. The proposed model may provide a theoretical framework for 
developing the refined model of the global pair-starved pulsar magnetosphere. 

Subject headings: theory — pulsars: general — stars: neutron 


1. INTRODUCTION 

It is believed that rotation-powered pulsars possess a magnetosphere in which charged 
particles are accelerated to relativistic Lorentz factors and generate a broad spectrum of 
pulsed emission (from radio, to IR, optical, X-ray and 7 -ray). The detailed physics of the 
global magnetosphere depends on various factors such as the efficiency of primary-particle 
(electron) ejection and acceleration near the neutron star (NS) surface, conditions for the 
occurrence and intensity of electron-positron pair creation, the spatial distribution of elec- 
tromagnetic fields and currents, and particle energy loss mechanisms. Although the vacuum 
magnetosphere model of Deutsch (1955) is, at present, the only available non-axisymmetric 
closed analytic model for the electromagnetic field of a rotating magnetic dipole, it is not 
an appropriate physical model of an active pulsar magnetosphere filled with charges and 
currents. Historically, the development of more realistic models began with the idealized 
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equation governing the structure of an axisymmetric pulsar magnetosphere (see e.g. Mestel 
1973; Scharlemann & Wagoner 1973; Michel 1973; Okamoto 1974; Mestel et al. 1979 and 
references therein). This equation, sometimes referred to as the pulsar equation, is the astro- 
physical counterpart of the force-free Grad-Shafranov equation (see the original publications 
by Grad [1967] & Shafranov [1966]). Further development was undertaken by Mestel (1999), 
Goodwin et al. (2004), Beskin et al. (1983), and Beskin, Kuznetsova & Rafikov (1998), 
and significant progress has been achieved in the numerical solution (see Contopoulos et 
al. 1999; Mestel 1999; Spitkovsky 2006; Timokhin 2006, 2007) of the pulsar equation in the 
ideal-MHD (E • B = 0) and force- free (neglect of particle inertia and pressure) approximation 
(see e.g. Arons 2004 for a brief theoretical overview). Although the force- free magnetosphere 
is probably a closer approximation to a real pulsar than the vacuum solution, it is still not a 
truly self-consistent model since the production of pair plasma requires particle acceleration 
and thus a break down of force-free conditions in some regions of the magnetosphere. Fur- 
thermore, it will apply to the younger pulsars that can readily supply the requisite charge 
through copious electron-positron pair cascades. 

A different approach to the theory of pulsar magnetospheres has been the simulations 
begun by Krause-Polstorff & Michel (1985, see also Smith et al. 2001, Petri et al. 2002), 
who studied how an axisymmetric rotating NS surrounded by vacuum fills with charge. 
They found that the evolution of a rotating, conducting sphere results in separated domes 
of charge over the magnetic poles and a torii of opposite charge in the equator, with vacuum 
between. Such a configuration is devoid of currents and charge outflow and is thus a dead 
pulsar. However, more recent 3D plasma simulations (Spitkovsky & Arons 2002; Biltzinger 
& Thielheim 2004; and Spitkovsky 2004, 2006, 2008) found that the charge domes are subject 
to diocotron instability, allowing charge from disrupted domes to fill the vacuum regions, 
opening up the possibility for a charge-filled magnetosphere. 

Studies of polar cap acceleration and pair cascades (Harding & Muslimov 2001, 2002, 
Harding, Muslimov & Zhang 2002) have found that only the youngest third of the known 
pulsar population (those that can produce pairs through curvature radiation) is capable of 
producing enough charges in pair cascades to completely screen the parallel electric field. 
The bulk of pulsars (those that produce pairs only through inverse-Compton radiation) can- 
not supply enough charge to screen the E\\ and will be “pair-starved” . The magnetospheres 
of pair-starved pulsars will thus lie in a very different regime from that of the force-free mag- 
netospheres. In our previous study (see Muslimov & Harding 2005; MHOS) we discussed the 
analytic solution pertaining to the pair-starved regime (Muslimov & Harding 2004; MH04) of 
particle outflow in the innermost magnetosphere of a NS with arbitrary pulsar obliquity, and 
derived the explicit expressions for the magnetic field and corresponding electric currents 
We assumed that the current of primary relativistic electrons was the only source of pertur- 
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bation for the originally pure dipole magnetic configuration. Our model implied therefore 
the occurrence of an accelerating electric field, E|| ^ 0, on open field lines. On the contrary, 
for the force-free and MHD models the electric current on open field lines is assumed to be 
very close to the Goldreich- Julian (GJ) current to ensure that E • B = 0. In this paper we 
propose a slightly different approach that lies between the MH05 and the force- free/MHD 
approximations described above, since it implies a self-consistent electric current for open 
field lines that is greater than the primary current but still smaller than the local GJ current, 
so that there is non-vanishing accelerating electric field along the open field lines. In this 
study we explore a fully analytic steady-state solution illustrating that the initially dipolar 
magnetic field undergoes growing topological change at increasingly higher altitudes by the 
self-consistent electric current that generates the 3D-monopole-type magnetic co nfi guration 
(the corresponding magnetic field strength, B, scales as 1/r 2 ). In our solution both the 
electric current and magnetic field are determined self-consistently, and it is assumed that 
near the polar cap surface the current matches the self-limited current calculated for the 
magnetic field determined by a pure dipole with small corrections from the 3D monopole. 

The paper is organized as follows. In § 2 we present a set of basic equations for the 
electromagnetic field that will be employed in our study. In § 3 we present our analytic 
solution for the magnetic field and currents. In § 4 we discuss the results of our study and 
summarize our main conclusions. 


2. Basic Equations 


The very general equations describing the electromagnetic field of a rotating NS in the 
Lab (inertial) frame are the first couple, 


V • B = 0, 

„ ^ 1 <9B 

V xE — — — , 
c ot 

and the second couple of Maxwell’s equations, 

V • E = 47 rp, 

„ „ IdE 4? r 

VxB = -- + -j, 
c at c 

where p and j are the electric charge and current densities, respectively; and all physical 
quantities, such as B, E, j, and p are defined in the Lab frame. 


(1) 

( 2 ) 

( 3 ) 

( 3 ) 


We will be searching for the steady-state solution to equations (l)-(4) in the domain of 
open-magnetic-field lines. In the steady state that is assumed to exist, the time derivatives in 
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equations (2) and (4) are determined by the proper set of special-relativistic transformations 
of the electromagnetic quantities in the case of rotational motion of open-field lines relative 
to the Lab frame. Moreover, it is important that in our model, even within the light cylinder, 
the open-magnetic-field lines may rotate slightly differentially without being wound up (see 
§3 below). To get the system of steady-state Maxwell’s equations, we can use the following 
special-relativistic transformations (see MH05) of partial time derivatives between the Lab 
frame (subscript “Lab”) and the frame of reference corotating with the open-field lines 
(subscript “corot”): 

'5B1 

OJc 


corot 


~dt) 


V x (u rot x B), 


Ltifc) 


dE 1 
dt J 


corot 


dE) 

dt J 


V x (u rot x E) + u rot V ■ E, 


( 5 ) 

(6) 


Lab 


dp\ 
dt J 


corot 


dp\ 
dt J 


+ u rot • V p, 


( 7 ) 


Lab 


where u rot (= ft x r, and 0 is the angular velocity that can, generally, be differential) is the 
linear rotational velocity of the open-field lines. 

We assume that in a steady state, the time derivatives in the LHS of equations (5)-(7) 
vanish (i.e. the electromagnetic fields “seen” by the observer corotating with the magnetic 
flux tube are stationary), so that Maxwell equations (2) and (4) can be rewritten in the 
following form 

V x E = -V x (/3 rot x B), 


47T 

V x B = V x (/3 rot x E) - /3 rot V • E + — j. 


( 8 ) 

( 9 ) 


where f3 rot = u rot /c. 

Also, the charge continuity equation, 


dp 

dt 


+ V-j = 0, 


with the help of relationship (7) takes the form 

Wot • Vp- V • j = 0. 

By combining equations (3) and (9), we get 

Urot 


( 10 ) 


:n( 


V x B 


4tt 

x E) = (j - pUrot)- 

c 


c 


( 12 ) 
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Thus, the steady-state solution for the domain of magnetosphere with open-field lines is 
determined by equations (1), (8), (11) and (12). Note that, since V-u rot = 0 and V • (pu rot ) = 
u ro t • Vp, equation (11) translates into 

V ■ (j - pu rot ) = 0. (13) 

We can ignore the second term in the LHS of equation (12) that is of order of (u mt /c) 2 <c 1. 
To complete the formulation we need to specify the current density j. We assume that 

j = j (1) + jrot +j E * B > (14) 

where is the density of the self-consistent electric current supporting the magnetic field 
B (1) (see formula [18] below) , j rot = pu rot is the density of the electric current generated by 
the bulk rotational motion of charges, and j ExB is the density of the E x B-drift current. 
In our previous paper (see MH05) we discussed the situation where the current j (1 ) was 
solely determined by the current of primary electrons in the dipole magnetic field. In this 
paper the current density j (1) matches the current density of primary electrons in the dipole 
field only at very small altitudes and provides significant distortion of initially dipole field 
at higher altitudes within the LC. Hence, in the present study, we may justifiably ignore 
much smaller contributions to the current j produced by the second-order correction to the 
primary electron current and by the E x B current, respectively, (see formulae [60]-[62] and 
[76]- [78] of MH05) in the RHS of equation (14). Because of the linearity of equation (12) the 
corrections to B d calculated in MH05 can be added to the final formulae for B (see equations 
[43] - [45] below). 

As in our previous study (see MH05), we shall use the magnetic spherical polar coordi- 
nates (x = r/R\ c , 9, 0; where R[ c = c/Q) with the polar axis along the magnetic moment of 
a NS. The magnetic coordinates are appropriate for modeling the magnetic field structure 
in the vicinity of the NS (at i< 1) with any obliquity. However, the magnetic coordinates 
are inconvenient in describing the effects of rotation on the global magnetosphere, simply 
because the rotation can break the symmetry with respect to the magnetic axis at x ~ 1 and 
for arbitrary obliquity. Intuitively, one may expect that, at higher altitudes, the magnetic 
coordinates get transformed into the spherical coordinates with the axis along the NS 
rotation axis that determines the global symmetry. Indeed, we can observe that (see e.g. 
formulae [73]- [75] in MH05), in magnetic coordinates, the effects of rotation (besides making 
the corresponding terms proportional to Q and having certain radial dependence) enter for- 
mulae for the magnetic field structure in an elegant way, through a single function, cosine 
of the angle between the NS rotation axis and radius-vector of a given point (in terms of 
magnetic polar coordinates), 


s = cos x cos 0 + sin x sin 9 cos 0, 


( 15 ) 
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and its derivatives over 9 and 0. Here, according to our notation, x i s the NS obliquity- 
angle. Note that s = 1 determines the symmetry axis in magnetic coordinates which is the 
NS rotation axis. In addition, the function s is a “toroidal potential” for the rotational linear 
velocity, /3 rot = u rot /c = (fir 2 /c) (Vs x n), where n is the unit vector in radial direction. 
We therefore choose to express our solution for which is in magnetic spherical polar 
coordinates, in terms of s. 


In our calculation the following couple of useful formulae will be employed, 

/ ds\ 2 1 / ds \ 2 _ 2 

\de) \Mj ~ ~ s • 


(16) 


and 


1 d 


ds 


1 d 2 s 


_ L/,5 \ X L/ o . . 

• a no \ sin ^7u) ) "• 2 7, Tua = ~2s. (17) 

sin 9 o9 \ o9 J sin 9 dtp 2 

These formulae illustrate that the function s is perfectly suitable for describing the effects 

of rotation (rotation-induced symmetry) in magnetic coordinates. 


For the sake of simplicity, we shall ignore all general-relativistic corrections (including the 
important effect of frame dragging ) throughout this paper. Since we are not going to discuss 
the polar cap electrodynamics, all these corrections may only unnecessarily complicate the 
resulting formulae. 


We assume, that the magnetic field in the innermost magnetosphere (well within the LC) 
can be presented as 


B = B (d) + eB (1) , 

(18) 

where 


Dd / 1 \ 

B d = —■ (cos 9e r + - sin 9eo ) 

(19) 

is a pure dipole magnetic held anchored into the NS and 


B (1) = B w (v,s,9,<t>) 

(20) 


is the magnetic held generated by the self-consistent electric currents in the domain of the 

magnetosphere with open held lines. The parameter e < 1 determines the strength of the B 1 '- 1 - 1 - 
component relative to the dipole one and also constrains the limiting value of the density of 
primary (electron) current at the polar cap surface (see condition [65] below). 

The Bd ( -component and the corresponding density of the self-consistent electric cur- 
rent are determined by the Maxwell’s equations (that, within the notations, coincide with 
equations [24j-[26] of MH05) 

€ V x B* 1 ’ = — j (1 >, 


( 21 ) 
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V • B (1 > = 0, 
V-j (1) = 0. 


( 22 ) 

(23) 


Equations (21)-(23) are the standard magnetostatic equations. Before we proceed with the 
analytic solution to these equations, we should specify the fundamental requirements to the 
expected solution. First, we will be solving equations (21)-(23) in magnetic spherical polar 
coordinates (r, 9,(p). Second, in our previous study (MH05) we have demonstrated that the 
generic current of primary electrons relativistically flowing along the dipole magnetic field 
of a NS is capable of generating the perturbation to the dipole magnetic field that has a 
l/r 2 -radial dependence. The fact that this current of primary electrons generates the 1/r 2 
correction to the magnetic field and matches the self-consistent current we derive in this 
paper at low altitudes (and in a small-angle approximation) is very suggestive. Also, the 
leading term in the charge density (all the way from the polar cap surface to the LC) scales 
as 1/r 3 which means that from basic dimensionality analysis Ffd) should scale as rp ~ 1/r 2 . 
So, we will be searching for the solution to equations (21)- (23) having the same radial 
dependence. Third, both the magnetic field and density of electric current are expected to 
depend on angular variables through the function s and its derivatives over 9 and (ft. This 
is important because in magnetic coordinates it is the function s that can be used in an 
elegant way to describe the effect of rotation and symmetry with respect to the rotation 
axis, and hence to ensure that the corresponding solution to Maxwell’s equations (21) - (23) 
will posses rotational symmetry. Finally, the ^-dependence will enter the solution (see also 
MH05) through the dimensionless radial coordinate, x = r/R\ c = flr/c. Physically, this 
means that we choose to express the solution for Bd) as (Dr/c) (Bp / t 7 3 )X, where X is some 
vector function (to be solved) that depends only on the function s and its derivatives over 9 
and (j). Likewise, the current density in equation (21) becomes automatically proportional to 
the flBo/27r?7 3 , which is just the amplitude of the current density of primary electrons. This 
general reasoning is sufficient to construct the exact solution for Bd) and jd) (note that jd) 
is just one of the three terms contributing to j, as discussed right after equation [14]) having 
the well-specified bulk properties. Thus, what remains to be determined is the exact form 
of angular dependence for our solution. 


3. Analytic Solution 

The general solution for B l lj satisfying solenoidality condition (22) and scaling as 1/j 
(and therefore can be dubbed as the 3D- mo nopole- type solution) reads 

B m _ BoZ-M 


4 


Vf c 


X * 


(24) 


(25) 


(i) _ 5 o 1 


B e = J 




1 (9s 

sin 9 dd) 


?(i) __ 


Bn 1 


T)f c X 2 


-o 9 P ( S ) 3 


1 9s 


sin# 99 


' ,,(s) I) ■ 


» P (S) = — s 2 ( 2? ) 

is uniquely determined by the solenoidality of Although this is a rigorous result, the 

singularity on the symmetry axis (s = 1) should be excluded from the solution domain for 
the simple physical reason that the symmetry axis is a current-free area (see also discussion 
following equation [48] below). In equations (24)-(26) x = r}/rji c , rj = r/R and r^ c = R\ c /R. 

Similarly, we can write the general solution for jW satisfying (23) as 

A (1) = 3 P(s), (28). 

dj 

® 1 + q,{s) ’ <29) 


;(1) _ *0 


p , , 1 Os t . . 0s\ 

q sin9 90 ~ q ^ 09 ) ’ 


where 


V (s) = y s [(1 - s 2 ) q p (s )] , 

as a direct consequence of the solenoidality of current (see [23]). 

By comparing the LHS of equation (21) with formulae (28)-(30), we find that 

n Bq 

Iq € 3 ? 

^vi 


q p = h\ 


Now we can present j (1 ) as a superposition of poloidal (superscript “p”) and toroidal (super- 
script “t”) components, 


and write 


j (,) =j p +j t 


j p = a p BP, 


j l = a‘ 
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where B p and B* are the poloidal and toroidal components of B (1) , and a p , a* are some 
scalar functions that will be determined below. 


By substituting expressions (32) - (34) into formulae (28) - (30) and making use of equations 
(21) and (35) - (37), we get 

(38) 


F = 

q p = h l = 


1 -s 2 ’ 

1 

(1 -s 2 ) 3 / 2 ’ 

2 


(39) 

(40) 

(41) 

2 7T (E \/l — S 2 

Thus, the components of the magnetic field (18) and electric current density (28)-(30) can 
be presented as 


cr 


(1-s 2 ) 2 ’ 

e 

4 7T X \/l — s 2 ’ 

e ft 


a 


(' 


cos 9 + e x 


1 — s 2 


B -®» 


1 sin 0 H — ——j 

2 1 — s 2 \9# Vl — s 2 sin# 90 


Bq e x 


1 9s 


1 


9s 


and 


rf 1 - s 2 V.sin 9 90 y/l - s 2 d9 J ’ 

, 0 ) = 


Jr 


ft Eg 


4l IT/ 3 (1 — S 2 ) 3 / 2 ’ 


(43) 

(44) 

(45) 

(46) 


? .(i) _ ftB$ 1 ^ 9s 2 1 9s ^ 

' Jf) 47 tt] 3 (1 - s 2 ) 3 / 2 \90 y'l - s 2 sin# 90/’ 

M 1 fj_ da 2 9s\ 

47T77 3 (1 — S 2 ) 3 / 2 \sin# 90 yT — 5 2 0# / ' 


(47) 

(48) 


One can see that formulae (43)-(48) have a singularity at the rotation axis, s = 1. The 
physical meaning of this singularity can be better understood by evaluating the GJ charge 
density (see formula [58] and discussion following formula [59] below) at the polar cap surface: 
the surface value of the GJ charge density decreases very near the rotation axis. This means 
that, the maximum charge density of electrons that can be ejected from the polar cap surface 
should also decrease toward the rotation axis. It is important that, on the magnetic field 
lines surrounding the rotation axis, the electrostatic potential is constant, and therefore the 
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accelerating electric field should vanish on these field lines. Thus, according to our solution, 
the polar cap region near the rotation axis is a current-free region that is also unfavorable 
for particle acceleration. Apparently, our analytic approach is not applicable to this region, 
and the field lines whose footpoints have s — s 0 = 1 should be excluded from consideration. 
We should also point out that, although our solution still implies some enhancement of the 
magnetic energy density towards the rotation axis (collimation of the field lines towards 
the rotation axis), the magnetic flux through the surface normal to the open field lines is 
constant and is by no means singular even at s — > 1. 

Now we should add to the (poloidal) current density of primary electrons, 

jpri m = ~ A ° ( cos 0e * + \ sin 6e <^ > ( 49 ) 


where 


Ao ~ cos x + x# sin x cos </>, 
z 


(50) 


in a small-angle approximation ( 8 1) that is valid at small altitudes above the pulsar 

polar cap surface. Note that, in MH05 the solution for B 1 ' 1 ) is determined by the density of 
primary electron current (equation [49]) only. In this paper we find that, in the pair-starved 
regime, there is a new type of the self-consistent solution that implies a gradual build-up of 
additional current (on the top of primary electron current) at increasingly higher altitudes 
while preserving the condition for the occurrence of non- vanishing accelerating electric field, 
E|| < 0. The magnitude of B^ given by formulae (43)-(45) is significantly larger than that 
of the corresponding component calculated in MH05. 


The total density of poloidal current can now be written as 

+ i p 

Jpnm 1 J ’ 


j t P ot=jn ^ iP 


(51) 


Note also that formulae (41), (42) imply that in the region (still within the LC) where B^ 
dominates over B !/t - , the magnetic surfaces (that do not intersect) are described by equation 

x y/l — s 2 = const. (52) 


On each magnetic surface, the angular velocity, Q, should satisfy the following condition 
(also known as the Ferraro’s isorotation law), 


B • Vfi = B p • Vfi = 0, 


(53) 


meaning that O is a function of w = x yl — s ' 2 only (in other words, the angular velocity 
is constant along each poloidal field line), thus ensuring that there is no winding-up of the 
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magnetic field lines. One can easily verify this by direct substitution. Formula (53) refers 
to the component B 1 ' 1 ^ or to the domain where dominates over B d . In the region 
where the magnetic field is mostly dipolar, the angular velocity of open field lines has two 
choices: (1) be constant everywhere (solid-body rotation); and (2) be constant along the 
field lines (in this case fl will be a function of sin 2 9/x that is a counterpart of function w for 
a pure dipole field). However, there is a problem with the second choice, simply because the 
dipole magnetic field (except the aligned case, x = 0) is not symmetric with respect to the 
rotation axis. Therefore, any differential rotation of the dipole field lines without crossing 
and entanglement would be topologically and physically prohibited, especially within the 
light cylinder. In this study we assume that, in the domain where the dipole field dominates, 
the angular velocity of open field lines is constant, 0 = = const (solid-body rotation), 

but as B^ begins to dominate Q — const is no longer required. 

To illustrate the effect of differential rotation of open field lines (in the domain where B ( b 
prevails), we may approximate Q(w) by the simple formula 

A = , (54) 

1 + aw 

where O 0 is the angular velocity of the NS, and a is a parameter (~ 1) to be determined by 
matching the solution for fi beyond the LC with that within the LC. This formula illustrates 
the simplest possible way for the angular velocity of open field lines to transit from the 
solid-body to differential rotation. It is, by no means, the actual solution for Q and will be 
used below just to examine and clarify some of the properties of our solution. 

Because our solution implicitly incorporates differential rotation of surfaces of constant mag- 
netic flux, we should be able to smoothly match the differential angular velocity beyond the 
LC with the nearly constant angular velocity well inside the LC (see e.g. the corresponding 
limits of formula [54]). It is important to note that in the case of differential rotation given 
e.g. by equation (54), our formalism can be extended even beyond the LC. More specifically, 
our formulae (5) - (7) are applicable for subluminal rotational linear velocity, Vtr \/ 1 — s 2 < c. 
By using equation (54) we can rewrite this fundamental constraint as w/(l + aw) < 1 that 
is actually valid both within the LC (w< 1) and beyond the LC (w > 1). This means that 
the problem of superluminal rotation gets automatically fixed: the magnetic field geometry 
changes to allow the differential rotation of field lines and, at the same time, ensure the 
subluminal rotation beyond the LC. 


For the region well within the LC (w « 1) , from formula (54) we can get, 

Q ~ fi 0 (l — aw). 


( 55 ) 


Now we can calculate the G.J charge density for the differentially rotating magnetic field 
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lines, 


Let us assume that 


Pgj = ■ (ftot X B )> 


n = F(w), (57) 

where F is some function of w (e.g. such as determined by formula [54]). 

By substituting (57) into formula (56) and performing necessary vector operations we 
can get 

n 15 In F\ QnB« fdF\ x 2 , . 


where 


Pgj Pgj ( ^ 


2 51n wj Attct ] 3 \dw J y/l — s 2 ’ 




is the canonical GJ charge density for a pure dipole field (in flat space-time), and 


A (9, 4>) — - [cos x(3 cos 2 9 — 1) + 3 sin x sin 9 cos 9 cos 


Note that the second term in equation (58) has the opposite sign (the differential velocity is 
likely to be attributed to the decrease in angular velocity of open field lines compared to the 
case of a solid-body rotation, i.e. dF/dw < 0) to the first term and is of order of ex 2 / y/l — s 2 . 
Consider the region around the rotation axis and the field lines whose footpoints have the 
values of so in the vicinity of so — 1- Given formula (55), that is perfectly applicable to the 
situation under discussion (i.e. where the field line geometry is essentially determined by the 
component and the differential rotation of open field lines becomes possible), and the 
fact that w should remain constant along the differentially-rotating poloidal field lines (see 
formula [52]), the second term in (58) should scale as (1 — s 2 )~ 3 / 2 , so that as s — s 0 — > 1 
it will tend to counterbalance the first term. The GJ charge density should significantly 
diminish towards the rotation axis before it changes sign near the very axis. Furthermore, 
the potential drop along those field lines tends to vanish since the potential is a f un ction of 
w, which is analogous to equation (53). This means that the ejection and acceleration of 
electrons from the small area of the polar cap surface surrounding the rotation axis shuts 
off, implying that our solution which assumes relativistic electron outflow is not applicable 
to this area containing the singularity. This is rather instructive: our solution formally 
contains a singularity but, at the same time, it indicates that near the singularity the initial 
model assumption (that electrons can be ejected from any point of the polar cap and then 
accelerated to relativistic velocity) gets violated. Therefore, the solution itself constrains 
its applicability by forcing the exclusion of the singular region from the solution domain. 
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A non-singular solution could be, in principle, obtained by incorporating the dynamics of 
particle acceleration from the polar cap surface into the model and by taking into account 
the feedback from the global structure of fields and currents (including the return current). 
Obviously, this kind of approach would dramatically complicate the model and prevent an 
analytic treatment. 

To estimate the space charge density, p, we should exploit the E|| = 0 boundary condition 
at the polar cap surface, or the approximate condition (valid within the accuracy of ~ 9$, 
where 9q is the canonical angular radius of the polar cap) that 

P~P G] at 77 = 1 . (61) 


Above the polar cap surface the observer “sees” the following charge density along the 
magnetic flux tube 

P=- c ~ 1 Jtof (62) 

By using formulae (43)-(45), (51), (61) and (62) the space charge density of relativistic 
electrons can be presented as 


Bj 

2ncri 3 rj 3 B 


A ° 2(1 — s 2 ) 3 / 2 ( A ^ ^Ao'v/l 1 _ g 2 j 


(63) 


Here A 0 and A are given by formula (50) and (60), respectively; and 


1 < 

2 (1 — Sq ) 3 / 2 



where s 0 ~ cos \ + >>o( sin \ cos d, ( = //(lj/ft, (0 < £ < 1), 0(1) is the magnetic colatitude 
of the footpoint of an open magnetic field line. 


Note that to constrain the model parameter e we should impose the condition 

« « (i - 4)2£, (M) 

implying that, at small altitudes above the polar cap surface, the corrections to p and 
p G] that are oc e are smaller than ~ 9% and do not affect the acceleration of primaries. 
Apparently, s 0 = 1 (rotation axis) is a geometrical singularity (see discussion following 
equation [48]), and the field lines with s 0 = 1 should be excluded from our solution. However, 
above the NS surface, in the situation where s o 1 the field lines get asymptotically (as 
r oo) collimated along the rotation axis (see equation [52]). The parameter e can be more 
accurately constrained by the condition Ej| = 0 at the surface. 

By using equations (43) - (45) one can “visualize” the field lines by taking cuts or projec- 
tions in different coordinate planes. However, we must admit that, without having a global 
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solution, one cannot reliably define the coordinates of footpoints of the last open field lines 
(we assume that the last open field lines are determined by the boundary of the pulsar dead 
zone where E-B = 0 provided that it gets formed). In this paper, for the sake of illustration, 
we assume that the last open field lines are emanating just from the rim of the canonical 
polar cap of radius r pc = 9$R which is, strictly speaking, valid only for the aligned case. 

In Figure 1 we presented a projection of the magnetic field lines in the X-Y plane as 
viewed down the pole for aligned rotator. The footpoints of the field lines have the same 
magnetic colatitude (£ = 0.35 and 0.4 for Figure la and lb, respectively) and magnetic 
azimuthal angles 0 = 15°, 30°, ..., 360°. For smaller values of £ there will be much stronger 
collimation of field lines. However, we caution that the region with very small values of £ 
may be free of particles and can be, therefore, excluded from our solution. Note that the 
unloaded (dipole) magnetic field lines emanating from this region can get deflected towards 
the rotation axis by ‘hoop stress’ generated by the surrounding loaded field lines. In Figures 
2-5 we depict the field line geometry in the meridional plane (cutting through the rotation 
and magnetic axes) for different values of obliquity angle, x = 0°, 30°, 60° and 90° and for 
the value of e — 0.01. To plot the field lines we used formulae (43)-(45) for the fixed values of 
magnetic azimuth, 0 = 0° and 0 = 180°. In Figures 2-5 we used the Cartesian X — Z plane 
where the positive Z-axis is along the magnetic axis, and the rotation axis (in the oblique 
cases) is pointing to the upper-right corner (in the orthogonal case the rotation axis is along 
the positive X-axis). The coordinates are scaled by the LC radius, R\ c . The magnetic 
colatitudes of the footpoints of all field lines range from 0.05 through 1.0 of the half-angular 
size of the canonical polar cap with the step of 0.05 (except the orthogonal case for which 
£ = 0.1, 0.2,. ..,1.0). From Figures 2-5 one can clearly see the effect of collimation of field 
lines along the rotation axis. The magnitude of the collimation depends on the parameter e 
and so (value of s for the field line footpoint): favorably-curved field lines tend to collimate 
much more than unfavorably-curved ones. For the aligned rotator (see Figure 2) the effect 
of collimation may be significant for most field lines emanating from the polar cap. Also, 
in the aligned case, the field lines experience substantial sweepback (see Figure la,b). It is 
important to note that, although our approximation is valid only within the canonical LC 
(up to ~ 0.5 — 0.8 of the LC radius), in Figures 2-6 we presented the magnetic field lines up 
to the cylindrical radius slightly greater than the LC radius (although some field lines close 
outside the LC, this is unlikely to happen in the global model, because here we have not 
included the E x B drift and rotation effect). In this way, we believe, the behavior of our 
solution could be better understood, and its deficiencies are more clearly exposed. Figures 
2-5 show that for small values of s the topology of magnetic field lines in the meridional 
plane (A' — Z plane) is very close to that of a dipole. 

Thus, the topology of the magnetic field lines in the proposed pair-starved model differs 
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from that of the force-free model (see most recent calculation by Kalapotharakos & Con- 
topoulos, 2008) in two major ways: the field lines in the pair-starved model are significantly 
more strongly collimated along the rotation axis, whereas the last closed field lines and adja- 
cent open lines may remain closed as opposed to the force-free model where these lines may 
form the“Y” shape structure due to centrifugal force and favor the formation of a current 
sheet. However, the detailed comparison of field geometries in the pair-starved regime with 
those in the force-free/MHD regime will be possible only after we present the closed global 
analytic solution. We expect that the combined effects of bulk rotation and E x B drift 
current (see formulae [60]-[62] and [76]-[78] in MH05) can significantly modify the topology 
of field lines in the (rotation) equatorial region. This and other issues will be addressed 
in a separate study aimed at the derivation of global analytic solution matching the inner 
solution discussed here. 


4. Discussion and Conclusions 

We have constructed a relatively simple self-consistent analytic solution for both the 
magnetic field and electric current in the inner magnetosphere (within the canonical LC) of 
a NS with arbitrary obliquity and for the pair-starved regime (see MH04). Our main assump- 
tion is that the so-called open magnetic field lines of a pulsar are loaded with relativistically 
outflowing charges. Also, in compliance with the pair-starved regime, we assume that in ad- 
dition to the primary electron current there is a supply of charges into the open-field region 
that is sufficient to feed the electric current capable of distorting the global pulsar magne- 
tosphere. In our model, the electric current (= — c|p|) is equal to the GJ current (= — c|p GJ |) 
only at the polar cap surface and varies with altitude (see formulae [51], [46]-[48]), remaining 
less than the local GJ current. Our solution is valid within the LC (say, for z\/l — s 2 <0.5) 
and implies the occurrence of non-vanishing accelerating electric field, Ey < 0. The self- 
consistent current, j (1 \ produces the effect of collimation of field lines along the rotation 
axis. For the oblique rotator, the effect of collimation is much more pronounced for the 
favorably curved (cos <fr > 0) than unfavorably curved (cos 6 < 0) field lines. The strongly 
collimated favorably-curved field lines (for which s — >■ 1) may not allow the continuous out- 
flow of electrons. Instead, these field lines may be loaded with electrons returning back to 
the stellar surface. It is also possible that near the rotation axis an acceleration-free zone can 
be established. For the unfavorably curved field lines the effect of collimation is much weaker 
but capable of moving the null-surface (defined by the condition ft ■ B = 0) to much higher 
altitudes or even completely eliminating it. At high altitudes and near the rotational equa- 
tor, where .i;vT - s 2 ~ 1, the E x B-drift current becomes important and should be taken 
into account. Qualitatively, this effect mitigates the collimation by causing the field lines to 
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flare further away from the rotation axis. Also, this effect plays an important role in current 
closure and should be incorporated into a truly global solution. In addition, our solution in- 
dicates that asymptotically, when becomes dominant, Ey — > 0, and therefore the regime 
of a steady-state acceleration may change to the force-free flow. Thus, our solution could be 
used to model the transition from the pair starved regime of the space-charge-limited flow 
to the force-free regime. 

The analytic solution we presented in this study can serve as a prototype for more 
refined global analytic models of the pair-starved regime. We expect that in a global solution 
the key unknown parameters, such as e.g. e will be more reliably constrained and related 
to the efficiency of particle acceleration. Such a global solution will explicitly include the 
second-order effects we neglected in this study. Hence, it would be interesting to explore, 
in a self-consistent manner, the effect of the bulk rotation and E x B-current and extend 
our analytic solution beyond the canonical light cylinder. Also, it would be interesting to 
investigate the effect of a possible change from the space-charge-limited to the force-free 
flow, as one goes from low to high altitudes, respectively (for the same pair-starved pulsar), 
and address the global current closure. In the next study, by using our analytic model, we 
will calculate the accelerating electric field within the domain of open magnetic field lines 
and update our previous calculations of Ey at very high altitudes, so that we will be able to 
explicitly take into account the effects of essentially non-dipolar global magnetic structure on 
the occurrence and geometry of high-altitude ‘slot gaps’ (Muslimov & Harding, 2003, 2004). 


We acknowledge support from the NASA Astrophysics Theory Program through the 
Universities Space Research Association. Also, we thank anonymous referee for the con- 
structive comments that helped to improve the manuscript. 
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Fig. 1 — The projections of the magnetic field lines onto the X — Y plane as viewed from 
the Z axis for the aligned case ( X = 0°). The footpoints of the field lines have magnetic 
colatitude of a) £ = 0.35 (left) and b) £ = 0.4 (right), and azimuths (p — 15°, 30°, ..., 360°. 
The parameter e — 0.01. 
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Fig. 5. The magnetic field lines in the meridional plane cutting through the rotation and 
magnetic axes for the orthogonal case (x = 90°). The footpoints of the field lines have 
magnetic colatitudes £ = 0.1, 0.2, ..., 1.0, and the parameter e = 0.01. 






